#!/usr/bin/env python
# -*- coding: utf-8 -*-

""" By Martin Senande-Rivera
    For Towards and atmosphere more favourable to firestorm development in Europe """

import os
import glob, sys
import numpy as np
import xarray as xr
import pandas as pd
from scipy import stats

GCM=sys.argv[1]  # Enter GCM
RCM=sys.argv[2]  # Enter RCM
rcp=sys.argv[3]  # Enter RCP scenario

path_outs='./'+RCM+'/'+GCM+'/'

# Number of days with FWI and Kindex above each threshold
FWI_Ki = xr.open_dataset(path_outs+'FWI-Ki_ndays_'+rcp+'.nc')['__xarray_dataarray_variable__']

# Linear trend function
def linear_trend(x):
    pf = np.polyfit(np.arange(1970.,2099.), x, 1)
    return xr.DataArray(pf[0])
# p-values function
def pvalue(x):
    r, p = stats.pearsonr(np.arange(1970.,2099.), x)
    return xr.DataArray(p)

# Linear trend of the number of days with FWI and Kindex above each threshold
trend = xr.apply_ufunc(linear_trend,  # function to apply
                         FWI_Ki,
                         input_core_dims=(('time',),),
                         vectorize=True) 
trend.to_netcdf(path_outs+'FWI-Ki_ndays_trend_'+rcp+'.nc')

# p-values of the number of days with FWI and Kindex above each threshold
pvalues = xr.apply_ufunc(pvalue,  # function to apply
                         FWI_Ki,
                         input_core_dims=(('time',),),
                         vectorize=True) 
pvalues.to_netcdf(path_outs+'FWI-Ki_ndays_pvalues_'+rcp+'.nc')
